﻿module Precession

open System
open Time
open Coordinates

let CalculateAxialPrecession (rightAscension:RightAscension, declination:Declination, currentEpoch:System.Int32, oldEpoch:System.Int32) =
    
    let degree = Convert.ToDouble(declination.Degree)
    let minute = Convert.ToDouble(declination.Minute)
    let second = Convert.ToDouble(declination.Second)

    let raH = Convert.ToDouble(rightAscension.Hour)
    let raM = Convert.ToDouble(rightAscension.Minute)
    let raS = Convert.ToDouble(rightAscension.Second)

    
    let Cdr = Math.PI / 180.0
    let Csr = Cdr / 3600.0;

    let Convh = Math.PI / 12.0
    let Convd = Math.PI / 180.0
    
    let decimalRightAscension = Coordinates.RightAscensionMinutesSecondsToDecimalRightAscension(raH, raM, raS)
    let decimalRightAscension = decimalRightAscension * Convh

    let decimalDeclination = Coordinates.DeclinationMinutesSecondsToDecimalDeclination(degree, minute, second)
    let decimalDeclination = decimalDeclination * Convd

    let a = Math.Cos(decimalDeclination);
    let array = Array2D.init 3 3 (fun x y -> 0.0)

    let x1 = Array.create 3 0.0
    x1.[0] <- a * Math.Cos(decimalRightAscension)
    x1.[1] <- a * Math.Sin(decimalRightAscension)
    x1.[2] <- Math.Sin(decimalDeclination)

    let deltaT = oldEpoch - currentEpoch
    let deltaT2 = Convert.ToDouble(deltaT)
    let t = 0.001 * deltaT2
    let deltaST = currentEpoch - 1900
    let deltaST2 = Convert.ToDouble(deltaST)
    let st = 0.001 * deltaST2
    let a = Csr * t * (23042.53 + st * (139.75 + 0.06 * st) + t * (30.23 - 0.27 * st + 18.0 * t))
    let b = Csr * t * t * (79.27 + 0.66 * st + 0.32 * t) + a
    let c = Csr * t * (20046.85 - st * (85.33 + 0.37 * st) + t * (-42.67 - 0.37 * st - 41.8 * t))
    let sina = Convert.ToDouble(Math.Sin(a))
    let sinb = Convert.ToDouble(Math.Sin(b))
    let sinc = Convert.ToDouble(Math.Sin(c))
    let cosa = Convert.ToDouble(Math.Cos(a))
    let cosb = Convert.ToDouble(Math.Cos(b))
    let cosc = Convert.ToDouble(Math.Cos(c))

    array.[0, 0] <- (cosa * cosb * cosc) - (sina * sinb);
    array.[0, 1] <- (-cosa * sinb) - (sina * cosb * cosc);
    array.[0, 2] <- -cosb * sinc;
    array.[1, 0] <- (sina * cosb) + (cosa * sinb * cosc);
    array.[1, 1] <- (cosa * cosb) - (sina * sinb * cosc);
    array.[1, 2] <- -sinb * sinc;
    array.[2, 0] <- cosa * sinc;
    array.[2, 1] <- -sina * sinc;
    array.[2, 2] <- cosc;

    let x2 = Array.create 3 0.0

    x2.[0] <- (array.[0, 0] * x1.[0]) + (array.[0, 1] * x1.[1]) + (array.[0, 2] * x1.[2])
    x2.[1] <- (array.[1, 0] * x1.[0]) + (array.[1, 1] * x1.[1]) + (array.[1, 2] * x1.[2])
    x2.[2] <- (array.[2, 0] * x1.[0]) + (array.[2, 1] * x1.[1]) + (array.[2, 2] * x1.[2])

    let mutable ra2 = Math.Atan2(x2.[1], x2.[0]);

    let conditional =
      if ra2 < 0.0 then ra2 <- (ra2 + (2.0 * Math.PI))
      else ra2 <- ra2

    let dec2 = Math.Asin(x2.[2])

    let ra2 = ra2/Convh
    let dec2 = dec2/Convd

    let coordinate = new Coordinate()

    coordinate.RightAscension <- Coordinates.DecimalRightAscensionToHoursMinutesSeconds(ra2)
    coordinate.Declination <- Coordinates.DecimalDegreesToDegreesMinutesSeconds(dec2)

    coordinate